#include "TMath.h"
#include <iostream>

using namespace std;


//=======================================================================================
void BinomialErrorInterval( Double_t Nsel = 8., Double_t Ntot = 10., Double_t * ErrorRegion = 0.){
//=======================================================================================
//---------------------------------------------------
// Purpose: Compute binomial error interval
// Author:  Ivo van Vulpen (Ivo.van.Vulpen@nikhef.nl)
// Date:    August 2007 
//---------------------------------------------------

  // Initialize return argument and check input arguments
  ErrorRegion[0] = 0.;
  ErrorRegion[1] = 0.;
  
  if ( Ntot <=0  || Nsel < 0 || Nsel > Ntot ){
    cout << "\n    ERROR: BinomialErrorInterval inputs:  Ntot <=0 || Nsel > Ntot \n" << endl;
    return;
  }

  // -------------------------
  // [1] Initialize parameters
  // -------------------------
  Double_t p_max          = Double_t(Nsel)/Double_t(Ntot); 
  Double_t OneSigmaRegion = TMath::Erf(1./TMath::Sqrt(2));
  Double_t p_right        = p_max;
  Double_t p_left         = p_max;
  Double_t frac_left      = 0.;
  Double_t frac_right     = 0.;

  // -------------------------------------------------------
  // [2] Compute error interval:
  //     a) Make approximation if necessary
  //     b) Integrate binomial to get 1 sigma error interval
  // -------------------------------------------------------

  //-- Test if we should make an approximation
  Double_t LogBetaFunction   = TMath::LnGamma(Nsel)+TMath::LnGamma((Ntot-Nsel))-TMath::LnGamma(Ntot);
  Int_t    MakeApproximation = ( LogBetaFunction < -700 ) ? 1 : 0; 

  //---------------------------------------------------------------
  //-- [2a] Make an approximation: error = sqrt( eff*(1-eff)/Ntot )
  //---------------------------------------------------------------
  if ( MakeApproximation ){
    Double_t error_estimate = TMath::Sqrt((p_max*(1.-p_max))/Double_t(Ntot));
    p_right    = p_max+error_estimate;  
    p_left     = p_max-error_estimate;  
  }

  //-------------------------------------------
  //-- [2b] No approximation: do the real thing
  //-------------------------------------------
  if ( !MakeApproximation ){

    Double_t frac_max    = TMath::BetaIncomplete(p_max, Nsel+1, (Ntot-Nsel)+1);
    //Int_t    Nsteps      = ( Ntot > 100 ) ? Ntot*1000 : Int_t(1e5);
    Int_t    Nsteps      = ( Ntot > 100 ) ? Ntot*100 : Int_t(1e4);
    Double_t dp          = 1./Double_t(Nsteps);  
    Double_t prob_left   = 1.;
    Double_t prob_right  = 1.;
    Int_t    KeepMoving  = 1;
    Int_t    Direction   = ( p_max < (1.-dp) ) ? 1 : -1; 

    while( KeepMoving ){
      
      //-- Move right
      if( Direction == 1 ){      
	p_right += dp;
	frac_right = TMath::BetaIncomplete(p_right, Nsel+1, (Ntot-Nsel)+1) - frac_max;
	prob_right = Nsel*TMath::Log(p_right) + (Ntot-Nsel)*TMath::Log(1.-p_right); 
      }
      
      //-- Move left
      if( Direction == -1 ){      
	p_left -= dp;
	frac_left = frac_max - TMath::BetaIncomplete(p_left, Nsel+1, (Ntot-Nsel)+1);
	prob_left = Nsel*TMath::Log(p_left) + (Ntot-Nsel)*TMath::Log(1.-p_left); 
      }
      
      //-- Decide which way to go    
      Direction = (prob_right < prob_left ) ? -1 : 1;
      if( Direction ==  1 && p_right > (1.-dp) ){ Direction = -1; } 
      if( Direction == -1 && p_left  < dp )     { Direction =  1; } 
      
      //-- Check if you can stop
      if( (frac_left+frac_right) > OneSigmaRegion ){
	KeepMoving = 0;
      }
      
    } // end KeepMoving

  }// end computing the full binomial error



  // ----------------------
  // [4] Return ErrorRegion
  // ----------------------
  ErrorRegion[0] = p_left;
  ErrorRegion[1] = p_right;

} // end BinomialErrorInterval()





//=============================================
void Example( Int_t Nsel = 8, Int_t Ntot = 10){
//=============================================

  Double_t ProbRegion[2];
  BinomialErrorInterval(Nsel,Ntot,ProbRegion);

  cout << " " << ProbRegion[0] << " "  << ProbRegion[1] << endl;

} // end Example()

